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Layered topological insulators, for example, Bi 2 Se 3 are optically hyperbolic materials in a range of 
THz frequencies. Such materials possess deeply subdiffractional, highly directional collective modes: 
hyperbolic phonon-polaritons. In thin crystals the dispersion of such modes is split into discrete 
subbands and is strongly influenced by electron surface states. If the surface states are doped, 
then hybrid collective modes result from coupling of the phonon-polaritons with surface plasmons. 
The strength of the hybridization can be controlled by an external gate that varies the chemical 
potential of the surface states. Momentum-dependence of the plasmon-phonon coupling leads to 
a polaritonic analog of the Goos-Hanchen effect. Directionality of the polaritonic rays and their 
tunable Goos-Hanchen shift are observable via THz nanoimaging. 


I. INTRODUCTION 

Bismuth-based topological insulators (TIs) have at¬ 
tracted much interest for their unusual electron surface 
states (SSs), which behave as massless Dirac fermions. 
However, bulk optical response of these compounds^^^® is 
also remarkable. The quintuple-layered structure of these 
materials causes a strong anisotropy of their phonon 
modes. The phonons that involve atomic displace¬ 
ments in the plane parallel to the basal plane (hence¬ 
forth, x-y or T-plane) have lower frequencies than ^ 2 ^, 
the c-axis (henceforth, z-axis) vibrations.^ For Bi 2 Se 3 , 
the dominant T- and z-axis phonon frequencies, 

= 64cm“^ =1.9 THz, 

= 135 cm-^ = 4.1 THz , 

differ more than twice. As a result, this and similar 
TIs can exhibit a giant anisotropy of the dielectric per¬ 
mittivity. There is a range of ui where the permittiv¬ 
ity tensor is indefinite: the real part of e^(w) is posi¬ 
tive, while that of e-*-(a;) is negative. Media with such 
characteristics are referred to as hyperbolic^®^^® because 
the isofrequency surfaces of their extraordinary rays in 
the momentum space k = are shaped as hy¬ 

perboloids [Fig. 1(a)]. In the THz domain, the widest 
band of frequencies where Bi 2 Se 3 behaves as a hyperbolic 
medium (HM) is between the aforementioned dominant 
frequencies, ^ < a; < however, other hyperbolic 

bands also exist in this TI (both at THz frequencies, see 
Sec. H, and at visible frequencies, see Ref. 19). It is im¬ 
portant that the approximate equation for the extraordi¬ 
nary isofrequency surfaces, 

(fc-)^ + (kv)^ {k^f _ 

e*(w) ^■'■(w) 

is valid up to jkj of the order of the inverse lattice con¬ 
stant. Accordingly, rays of momenta jkj greatly exceed¬ 
ing the free-space photon momentum w/c can propa¬ 
gate through hyperbolic materials without evanescent de¬ 
cay. At such k the hyperboloids can be further approx¬ 
imated by cones, which means that the group velocity 



FIG. 1. (Color online) (a) Hyperboloidal isofrequency sur¬ 
faces of HP^s for two frequencies u>i and L 02 { 0 J 2 > toi). The 
asymptote angle 6 with respect to the k^-k^ plane is shown; 
the group velocity v makes the same angle with respect to 
the fc^-axis. (b) Model geometry: a TI slab of thickness d 
sandwiched between a substrate of permittivity ts and a su¬ 
perstrata of permittivity cq. The two thin (orange) layers 
represent the top and the bottom surfaces states. 


V = duj/dk of the rays makes a fixed angle 0 (or —6) 
with respect to the z-axis, with 


ta,n9{oj) 


[e"(w)]i/2 ’ 


( 3 ) 


see Fig. 1(a). We refer to these deeply subdiffractional, 
highly directional modes as the hyperbolic phonon po¬ 
laritons (HPP or HP^, for short). 

Our interest to HP^ of TIs is stimulated by recent 
discovery^^’^^ and further exploration of similar collec¬ 
tive modes in other systems such as hexagonal boron 
nitride^^^^® (hBN) and hBN covered by graphene^®^^® 
(hBN/G). There is a close analogy between these sys¬ 
tems. In fact, except for the difference in the number 
of Dirac cones {N = 1 vs. N = A) and the frequency 
range where the hyperbolic response occurs (THz vs. 
mid-infrared), the electrodynamics of longitudinal collec¬ 
tive modes of Bi 2 Se 3 and hBN/G structures is qualita¬ 
tively the same. (The analogy is the most faithful when 
graphene and hBN are rotationally misaligned; other- 
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wise, their collective modes are modified by the moire 
superlattice effects. 

The main goal of this paper is to investigate the in¬ 
teraction of HP^ with the Dirac plasmons of the topo¬ 
logical SS. The latter dominate the charge (and cur¬ 
rent) density response of the system at frequencies out¬ 
side the hyperbolic band where HP^ are absent. Dirac 
plasmons have been extensively studied in previous lit- 
erature®d3d4, 30-44 both TI and graphene. The basic 
properties of the Dirac plasmons can be introduced on the 
example of a hypothetical TI material with a frequency- 
independent permittivity e'^ > 0 and the permittivity 
£■*■( 0 ;) dominated by a single phonon mode. Such an ide¬ 
alized material is hyperbolic in a single frequency interval 
ojto < u! < ujio where £"'■( 0 ;) < 0. Its Dirac plasmons ex¬ 
ist at w < LOto and uj > ujio where > 0. In the 

setup shown in Fig. 1(b), where the TI slab borders me¬ 
dia of constant permittivities £0 > 0 and £s > 0 , there are 
two plasmon modes. At large enough in-plane momenta 
q = [(A:“)^-I-(fc^)^]^/^ these modes are confined to the op¬ 
posite interfaces and electromagnetically decoupled. In 
the relevant range of momenta g < g*, the dispersion of 
the plasmon bound to the top interface is given by 

g(^) - ^ ’ hu}<^\^i\, (4) 


where 


£l(u;) = [£^(a;)]l/M£^(u;)]l/^ (5) 

is the effective permittivity of the TI and /r is the chem¬ 
ical potential of the SSs measured from the Dirac point. 
At frequencies far below coto or far above uiio, function 
£ 1 ( 0 ;) can be approximated by a real constant, which 
yields uj oc ^/q. This typical two-dimensional (2D) plas¬ 
mon dispersion describes the low-frequency part of the 
full curve sketched in Fig. 2(a). The plasmon dispersion 
for the bottom interface is obtained by replacing £9 with 
£s (unless £s ^ £ 0 , in which case the range g > g* is rel¬ 
evant where the dispersion is approximately linear, see 
Sec. IIIB). 

Equation (4) implies that the nature of the plasmon 
modes should change drastically when ui enters the hy¬ 
perbolic frequency band where £ 1 ( 0 ;) [Eq. (5)] is imagi¬ 
nary and strongly w-dependent. This equation predicts 
a complex g, which suggests that the Dirac plasmons be¬ 
come leaky modes that rapidly decay into the HP^ bulk 
continuum. However, this is not quite correct. We will 
show that nonleaky, i.e., propagating modes can survive 
in thin enough TI slabs where the HP^ continuum is bro¬ 
ken into discrete subbands of waveguide modes. The lat¬ 
ter hybridize with plasmons to form hyperbolic plasmon 
phonon polaritons (HPPP or HP^, for short), the pri¬ 
mary target of our investigation, see Pigs. 2(b) and (c). 
We explore the following properties and manifestations 
of the collective charge modes of the TIs: i) the mode 
dispersion in the momentum-frequency space, ii) the de¬ 
pendence of such dispersions on the surface doping and 


the thickness of the slab, iii) the unusual real-space dy¬ 
namics of the HP^ rays, including a polaritonic analog of 
the Goos-Hanchen (GH) effect. 

The remainder of the paper is organized as follows. 
In Sec. II we specify the model and the basic equations. 
In Sec. Ill we present our results for the dispersion of 
the three different types of collective modes (plasmons, 
HP^s, and HP^s). In Sec. IV, which is the centerpiece 
of this work, we discuss waveguiding and launching of 
the HP^ modes and also their tunable GH shifts. We 
explain how these phenomena can be probed experimen- 



FIG. 2. (Color online) Schematic illustrations of the collective 
mode spectra in idealized model systems, (a) The plasmon 
dispersion of Dirac fermions confined to the interface of two 
bulk media of constant positive permittivity eo and e^. The 
dispersion crosses over from ui ~ to ca ~ ng at a 

characteristic momentum g, [Eq. (26)]. The shaded areas 
indicate the electron-hole continua where the plasmons (and 
any other charged collective modes) are damped, (b) The 
dispersion of hybrid HP^ modes for a slab of a hypothetical 
TI material that has a single in-plane phonon mode at wto and 
constant > 0. Permittivity is negative at Wto < cu < ojjo 
and positive at other w. The dotted boundary corresponds to 
the dotted line in (a). Outside the band Wto < cu < wuo, only 
plasmonic modes 0 and 1 exist. In the degenerate case eo = 
Cs they correspond to the symmetric (s) and antisymmetric 
(a) combinations of the top and bottom interface plasmons. 
Inside that band, multiple branches of HP^ are formed due to 
hybridization of the plasmons with the HP^ waveguide modes. 
The frequencies of all the branches other than 0 and 1 tend 
to wio at large momenta, (c) Schematic in-plane electric field 
profiles of the first few HP^ modes (thick curves). The number 
of nodes in each profile (the points where they cross with the 
vertical lines = 0) is equal to the modal index. 
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tally using the imaging capabilities of the scattering-type 
scanning near-held optical microscopy (s-SNOM).'^^d8 
Sec. V we give concluding remarks and an outlook for 
the future. Finally, in Appendix we discuss signatures 
of the phonon-plasmon coupling measurable by the s- 
SNOM operating in the spectroscopic mode. 


II. MODEL 


Our model for the bulk permittivities of the TI is 

e“(w) = + XI "2 ' .Fa- ’ a = -L, z . (6) 


i=i.2 


- Z7“a; 


In the case of Bi2Se3, we choose the parameters based 

on available experimentapdd theoretical^ literature 

as follows: = 29, = 17.4, = 64cm“^, = 

704 cm“^, a;^2 = 125cm“^, ujp 2 = 55cm“^, i = 
135cm“^, Wpi = 283cm“^, 0Jtg2 — 154cm“^, a;p2 = 
156cm“^, and 7“ = 3.5cm“^. [Note that ^ and ^ 
were already listed in Eq. (1).] The real parts of functions 
£-*-(0;) and e^(w) are plotted in Fig. 3. The regions where 
at least one of them is negative are shaded. They include 
region A, 3 < uj < oJto^u where Bi2Se3 is a HM of type 
II (3fiee2 > 0,3fiee-*- < 0); region C, ^ < w < 163 cm“^ 
where it is a HM of type I (3fieez < 0, 3fiee-*- > 0), and 
region B, i < co < 146 cm“^, where it exhibits the 

Reststrahlen behavior i^ee^ < 0, Sftee-*- < 0). Since 
regions B and C are narrow, in our discussion of HP^ and 
HP^ modes we focus on region A. In this discussion we 
often refer to hBN as an example of a simpler material. 
The type II hyperbolic band of hBN is bounded by the 
frequencies^®’^^ 

ujto = 1376cm“^, w/o = 1614cm“^. (7) 

In this band of hBN can be modelled similar to 

Eq. (6) but using a single Lorentzian oscillator while 
can be considered w-independent and positive. 

In the case of Bi2Se3, we also have to specify our as¬ 
sumptions about the electronic response. We consider 
only frequencies smaller than the bulk gap 0.3 eV of 
Bi2Se3 at which the electronic contribution to the per¬ 
mittivities [included in Eq. (6) via is purely real. 
Additionally, we assume that the valence bulk band is 
completely filled, the conduction one is empty, with no 
free carriers present in the bulk. However, such carri¬ 
ers populate the gapless SS described by the massless 
2D Dirac equation. The chemical potential which is 
located inside the bulk band gap, determines the dop¬ 
ing of these SS. We ignore any virtual or real electronic 
transitions between the surface and the bulk states. 

The fundamental current/density response functions of 
the SS are the sheet conductivity a and polarizability P, 
which are related in the standard way: 

= ^e^P{q,oj). (8) 

r 



FIG. 3. (Color online) The real parts of the tangential and 
axial permittivities of Bi 2 Se 3 . The sign changes of the per¬ 
mittivities are due to the Eu and A 2 U phonons. Surface- and 
bulk-confined collective modes exist inside the spectral regions 
where at least one of the permittivitties is negative. They in¬ 
clude type II hyperbolic region A (Kee^ < 0, IRee^ > 0), 
Reststrahlen region B (5Ree^,5Ree^ < 0), and type I hyper¬ 
bolic region C (IRee^ > 0, Kee^ < 0). 


Within the random-phase approximation for Dirac 
fermions, P{q,uj) can be computed^®’^° analytically: 


= - 


Nkp 


iN 


2'Kfiv \QttTlv ^q^ — 


G 


kijj 2kp 


-G 


k^ - 2kp 


q J \ q 

G{x) = ix\l 1 — — i arccos x . 


( 9 ) 


Here the branch cut for the square root and logarithm 
functions is the negative real semi-axis, k^^ is defined by 
kui = {to + tye)/^, phenomenological parameter 7e > 0 is 
the electron scattering rate, v is the Fermi velocity, and 
kp = \^\l{hv) is the Fermi momentum. Equation (9) is a 
good approximation at small At large doping, trigonal 
warping^^ and other details of realistic band-structure'*^ 
should be included. Since the above formula is a bit cum¬ 
bersome, it may be helpful to mention some properties 
of a{q,ui). For example, if 7e = -1-0, the real part of 
a{q,u}) is nonvanishing only inside the two shaded areas 
in Fig. 2(a), which together form the so-called electron- 
hole continuum.(This real part is a measure of dis¬ 
sipation, i.e., Landau damping.) For doped system at 
small momenta and frequencies, q, k^j kp, the expres¬ 
sion for the conductivity can be reduced to 


Ne^ kp ik,^ 

27r?i ,^q2 _ 1.2 ^q2 _ j.2 


( 10 ) 


At q <C uj/v, it further simplifies to the Drude formula 


Ne^ 1^1 

47r?i^ 7e - feu ’ 


(j{q,u}) 




( 11 ) 
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For an undoped system, one finds instead 


(j{q,uj) 


N iki^ 


N 

16 ¥ ’ 


LU 

V 


( 12 ) 

(13) 


In order to find the dispersion of the collective modes 
of the TI slab we use two computational methods. One 
method, which is advantageous for deriving analytical 
results, is to look for the poles of the response function 
rp(q,uj). This function is the total P- (also known as 
the TM-) polarization reflectivity of the system measured 
when an external held is incident from the medium la¬ 
beled “eo” in Fig. 1(b). It must be immediately clarihed 
that rp{q,uj) has no poles at simultaneously real q and 
w if the dissipation parameters 7 and 7e are nonzero. At 
least one of these arguments must be complex. When¬ 
ever one refers to the dispersion relation of a mode, one 
means the relation between the real parts of q and lo. The 
other method, which is especially convenient for numeri¬ 
cal simulations, is to identify the sought dispersion curves 
with the maxima of Amrp(q,a;) at real arguments. As 
long as the imaginary parts of q and lu (which give in¬ 
formation about the propagation length and lifetime of 
the mode) are small, both methods give the same disper¬ 
sions. An extra benefit of working with real q and w is 
that the corresponding rp{q,u!) is the input for further 
calculations we discuss in Appendix A where we model 
s-SNOM experiments for the system in hand. 

Our procedure for calculating function rp{q, uj) can be 
explained as follows. Taking a more general view for a 
moment, we regard the entire system including the sub¬ 
strate and superstrata as a stack of layers j = 0,1 ,..., M 
of thickness dj , tangential permittivity , and axial per¬ 
mittivity €j. (In the present case, M = 2, the TI slab is 
layer j = 1 and di = d.) Additionally, we assume that 
the interface of the layers j and j + l possesses the sheet 
conductivity Ujj+i. We observe that the P-polarization 
reflectivity of j,j + 1 interface in isolation is given 

by the formula (see, e.g.. Ref. 27) 




Qj+i Qj T 


Qj+i + Qj + 


47r 

77'^ij-i-i 
LU _ 

47r 

77'^jj-l-i 

CJ 


Qj — 




(14) 


(15) 


where fc| and q are, respectively, the axial and the tan¬ 
gential momenta inside layer j. Let rj be the reflectivity 
of a subsystem composed of layers j, ■ ■ ■, M. By this def¬ 
inition, rM-i = The crucial point is that the 

desired rp = ro can be found by the backward recursion 


Q-rj,j+i)Q-rj+i,j)rj+i 
rjpijrj+i - exp{-2ikj+idj+i) ’ 


(16) 


where rjpij is the right-hand hand of Eq. (14) with Qj 
and Qj+i interchanged. For M = 2, one recursion step 


suffices, which gives us, after some algebra, 

?'i2(?'oi + rio - 1) - roi exp(-2i/cidi) 

Tp = - 

riori 2 - exp{-2ikidi) 

Hence, function rp{q,ui) has poles whenever 


(17) 


rio[q, uj)ri 2 [q, uj) = exp {-2ikld) . (18) 

For large in-plane momenta q ^ (w/c) max we 

can use the approximations k^ q tan 6 and 


eo - ei 


rio - 


9top 


eo + ei - ^ 

ytop 


9top — 


IW 


27r(Ji 


top 


(19) 


where crtop = atop(9, w) is the sheet conductivity of the 
SS at the top interface. Let us also define the “phase 
shifts” (()top and (j)hot for inner reflections from the top 
and bottom interfaces, respectively: rio = — exp{2i(j)top), 
ri 2 = — exp(2i(()bot)- Note that in general (j)top and (fihot 
are complex numbers. Specifically, we take 


<(>top = arctan 



4>hot = arctan 


. Ss 
i — 

ei 


1-¥ j 

Co 9top / 
QhotJ 


( 20 ) 

( 21 ) 


where the standard dehnition of arctan z is assumed, 
with the branch cuts (—ioo,—z), (i,ioo) in the complex- 
z plane; 9bot is defined analogously to qtop but with the 
sheet conductivity CTbot of the bottom SS instead of atop. 
Equation (18) can now be transformed to 

2 

qn = + (j)top + 4>hot) , 6 = 2dtan6, (22) 

where the integer subscript n labels possible multi¬ 
ple solutions. Admissible n must satisfy the condition 
Amrp(9„,a;) > 0. Our numerical results for rp com¬ 
puted from Eq. (17) and analytic approximations for the 
solutions of Eq. (22) are presented in Sec. III. 


III. COLLECTIVE MODE DISPERSIONS 

The false color maps of function QTarp{q,uj) provide 
a convenient visualization of the collective mode spectra. 
Examples of such maps computed for Bi2Se3 slabs are 
presented in the bottom row of Fig. 4. Their counterparts 
for graphene-hBN-graphene (G/hBN/G) structures are 
shown in the top row to facilitate the interpretation. The 
bright lines in Fig. 4 are the dispersion curves of the 
collective modes. The apparent widths of those lines give 
an idea how damped the modes are. Below we discuss 
these results in more detail. 


A. Hyperbolic waveguide modes 

Figures 4(a) and 4(d) depict the Qmrp maps for, re¬ 
spectively, G/hBN/G and Bi2Se3 slabs, when they are 
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FIG. 4. (Color online) Collective mode dispersions of graphene-hBN-graphene (G/hBN/G) and Bi 2 Se 3 slabs rendered using 
the false color maps of Srmrp. The parameters of the calculation for G/hBN/G are: (a) d = 60 nm, p = 0, (b) d = 60 nm, 
= 0.29eV, (c) d — 30nm, p = 0.29eV. The other parameters are v = 1.00 x 10® cm/s, 7 e = 1.00THz, eo = 1, and es = 1.5. 
The parameters of the calculation for Bi 2 Se 3 are: (d) d — 120 nm, p = 0, (e) d = 120 nm, /r = 0.29 eV, (f) d = 60 nm, 
/i = 0.29 eV. In these three plots v = 0.623 x 10® cm/s, 7 e = 1 THz, eg = 1, and = 10. Equal doping of the top and bottom 
SS is assumed. The vertical dashed lines indicate a characteristic momentum probed by the s-SNOM experiments simulated 
in Fig. 7 below. 


undoped, /i = 0. No Dirac plasmons exist in such sys¬ 
tems, so that the collective modes are limited to HP^s. 
In Fig. 4(a) we see a single family of such modes whereas 
in 4(d) one can actually distinguish three of them. Let 
us start with the former, simpler case. The key to un¬ 
derstanding the nature of these modes is that inside the 
hyperbolic band LOto < OJ < loio the z-axis momentum 
kl ~ gtan0 of the modes is nearly real. Hence, the HP^s 
form standing waves inside the slab. The integer n in 
Eq. (22) corresponds to the number of nodes of these 
waves, see Fig. 2(c). For G/hBN/G the requisite condi¬ 
tion SJmrp > 0 is satisfied by all nonegative integers n 
due to the fact that 9m tan 0 > 0. This inequality also 
ensures that Qmq > 0. An analytical approximation for 
the dispersion curves of an undoped slab is obtained by 
neglecting the fractions q/qtop, q/qhot in Eqs. (20), (21), 
in which case Eq. (22) yields q{uj) directly. Within this 
approximation, momenta qn at given oj are equidistant: 

27r TT 1 

(23) 

The dispersion of the HP^ waveguide modes is dominated 
by the factor 1/tan0(a;) in Eqs. (22), (23), which, if all 
damping is neglected, changes from zero to inhnity as w 
increases from ujto to ivio- This is precisely what we see in 
Fig. 4(a): all the dispersion curves start at uJto at g = 0 
and increase toward ujio at large q. 


Equation (23) is general and it applies to Bi2Se3 as 
well. The three families of collective modes seen in 
Fig. 4(d), belong to the spectral regions A, B, and G of 
Fig. 3. In region A, which is the widest of the three, we 
see a set of HP^ modes very similar to those in Fig. 4(a). 
They start at uJto,i = 64cm“^ at g = 0 and monotonically 
increase toward uJto ,2 = 135 cm“^ at large g. In region 
C, 154 < u! (cm“^) < 163, we again hnd a family of HP^ 
modes but this time with a negative dispersion. This be¬ 
havior is typical of type I HM (3fiee-*- > 0, 3fiee^ < 0). 
The shape of the dispersion can be understood noticing 
that the real part of 1/tan0(a;) is positive, varying from 
oo to 0 (if the phonon damping 7“ is neglected) while ad¬ 
missible n are now n < 0. [In hBN, this type I behavior 
is also realized^^’^^’^^ but the corresponding frequency 
range is below the axis cutoff in Fig. 4(a).] Lastly, in 
region B, 135 < w (cm“^) < 146, function tan0(a;) is al¬ 
most purely imaginary, which implies that the collective 
modes do not form standing waves but are exponentially 
confined to the interfaces. Also, there are only two such 
modes, n = 0 and n = 1. In this respect these surface- 
bound HP^ modes are similar to the Dirac plasmons, see 
Sec. I above and Sec. IIIB below. However, their disper¬ 
sion is completely different from those of the plasmons, 
e.g., the dispersion of the upper (n = 1) mode has a neg¬ 
ative slope, see Fig. 4(d). Similar collective excitations 
have been studied in literature devoted to other systems. 
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e.g., anisotropic superconductors,®^ which can be con¬ 
sulted for details and references. Due to narrowness of 
regions B and C, some of the described features may be 
difficult to see in Fig. 4(a) and probably challenging to 
observe in experiments. For this reasons, we will mostly 
refrain from discussing regions B and C further. 

One implication of Eq. (23) is that the HP^ dispersion 
is widely tunable: the scaling law cx d~^ provides 
a practical way to engineer a desired wavelength of the 
waveguide modes simply by tailoring the slab thickness d, 
as has been previously demonstrated using hBN slabs. 


B. Surface plasmons 

Examples of the collective mode spectra at hnite dop¬ 
ing are shown in Eig. 4(b, c) for G/hBN/G and 4(e, f) 
for Bi 2 Se 3 . The spectra are dramatically different inside 
and outside the hyperbolic frequency bands. A key to un¬ 
derstanding this difference is again the value of the mo¬ 
mentum k\ ~ gtand(a;). Outside the hyperbolic bands, 
it is almost purely imaginary, and so the collective exci¬ 
tations are exponentially confined to the surfaces of the 
slab. These surface modes are the Dirac plasmons intro¬ 
duced in Sec. I. Having in mind applications to near-held 
experiments, we are particularly interested in momenta 
q of the order of a few times 10 ®cm“^, i.e., the region 
nearby the dashed lines q = 0.0025 nm“^ in Fig. 4. If ei 
is real, there are at most two solutions of Eq. (22), one 
for n = 0 and the other for n = 1. However, the distinct 
n = 1 dispersion curves are visible only in Eig. 4(b, c) 
for G/hBN/G and none of them is close enough to the 
range of q we are interested in. Therefore, we focus on 
the n = 0 branch. 

The shape of the plasmon dispersion curves in TI slabs 
and double-layer graphene systems was a subject of many 
previous theoretical studies®^’®^’^®’^®’^^ whose basic con¬ 
clusions are reproduced by the following analysis. To 
the right of the dashed lines in Fig. 4(b, e) and for 
d ~ 100 nm, the dimensionless product 2kfd = q5 is typ¬ 
ically large by absolute value and almost purely imagi¬ 
nary. This implies that the plasmons of the two interfaces 
are decoupled. Taking into account that eg < Cs and 
Qtop = Qbot in Fig. 4, one can show that the dispersion of 
the n = 0 mode is controlled by the properties of the top 
interface. In the hrst approximation this dispersion can 
be obtained setting (/top —t —ioo, which yields 

90 « 9top , qo > |(5r^ (24) 

For /r = 0, momentum (/top = 9 top((Zo,w) is imaginary, 
cf. Eqs. (13) and (19). Hence, for real ei, Eq. (24) has no 
real solutions: as already mentioned, undoped SS do not 
support plasmons. Indeed, Figs. 4(a) and (d) contain no 
bright lines outside the hyperbolic bands. On the other 
hand, if // 0, we can use Eq. (11) to transform Eq. (24) 

to Eq. (4), which predicts a parabolic dispersion curve 
to cx y/q if ei is constant. Such parabolas are seen in 


the upper halves of Figs. 4(b, c) and (e, f) although they 
appear rectilinear because of the restricted range of q. 

As smaller momenta Eq. (24) no longer holds. The 
correct approximation for the n = 0 mode is obtained by 
setting the left-hand side of Eq. (22) to zero. This yields 
c/top — 4^hot and 


Eq + Es 1 

2 9top + 9bot 


2 

N 


Co + Cs 




(25) 


Thus, both the \ow-q and high-(/ parts of the n = 0 dis¬ 
persion curve are parabolic but with different curvatures. 
The crossover between these two parabolas occurs via a 
rapid increase of £"'■( 0 ;), and so, ei(w) at frequencies im¬ 
mediately above the hyperbolic bands. It takes place at 
w > 1614cm“^ for G/hBN/G and w > 163cm“^ for 
Bi 2 Se 3 , which generates the inflection points seen on the 
curves in, respectively, Fig. 4(b, c) and (e, f). 

As indicated schematically in Fig. 2(a), at very large 
q the plasmon dispersion should have another inflec¬ 
tion point. Using the more accurate Eq. (10) instead 
of Eq. (11), we Hnd the following analytical result for the 
frequency of the n = 0 mode as a function of q: 


U}{q) 


q + q* 

\J\ -I- (2(/*/(/) 


2e2 Nkp 

hv Co + €i 


(26) 


This equation predicts a crossover from the parabolic to 
the linear dispersion u) ~vq above q = q*. However, this 
occurs far outside the plot range of Fig. 4. 

Returning to Eq. (25), we notice that it does not con¬ 
tain the bulk permittivities. Hence, it should continue 
to hold for a range of ui inside the hyperbolic bands. A 
physical picture of this mode [“0(s)” in Fig. 2(c)] is in- 
phase oscillations of the charges of both Dirac fermion 
layers, i.e., the system behaving as a single 2D layer with 
the combined oscillator strength. As oj decreases further 
into the hyperbolic bands, the length scale |(5| increases. 
The strength of the inequality (/o]/] «C 1 and so the ac¬ 
curacy of Eq. (25) becomes progressively lower [in fact, 
Eq. (27) below gives a better approximation]. At w = toto 
for G/hBN/G and similarly, at w = w// 3 for Bi 2 Se 3 , this 
inequality is violated completely, which is consistent with 
the termination of these branches at (/ = 0 in Figs. 4(b) 
and (e). Similar analysis can be applied to Figs. 4(c) 
and (f) where d is twice smaller than in, respectively. 
Figs. 4(b) and (e). Because of that, the plasmon disper¬ 
sion in the region (/]/] < 1 is shifted to smaller q. The 
dispersions in the large-(/ regions are virtually unaffected 
since the stronger surface confinement of the plasmons 
makes them insensitive to d. 

One qualitative difference between G/hBN/G and 
Bi 2 Se 3 is the richer phonon spectrum of the latter. This 
leads to the avoided crossings of the plasmon branch with 
the dispersion lines of the HP^ modes in regions B and 
C of Bi 2 Se 3 , cf. Fig. 4(b, c) and 4(e, f). The small shifts 
caused by those crossings are somewhat masked by the 
considerable linewidth of the n = 0 line due to rela¬ 
tively stronger phonon damping. In turn, higher elec¬ 
tronic damping rate 7 e ~ w/; 3 due to disorder scattering 
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in Bi 2 Se 3 effectively eliminates the plasmon excitations 
in the lower spectral region a; < see Fig. 4(e, f). 

Therefore, we do not discuss it here. 


C. Hybrid modes 


From now on we turn to the subject of our primary in¬ 
terest, the hyperbolic collective modes of a doped TI. 
In this short section we address their dispersion law. 
Comparing Fig. 4(d) for = 0 with Fig. 4(e, f) for 
/j, > 0, we observe significant shifts in the dispersion of 
the n = 0 mode in the upper half of the hyperbolic band 
I < uj < LOfg I of Bi 2 Se 3 . Similar shifts are seen in 
hBN near uJio, cf. Fig. 4(a) with Fig. 4(b, c). These 
shifts result from hybridization of HP^ and Dirac plas- 
mons into combined HP^ waveguide modes. In general, 
calculation of these shifts requires solving Eq. (22) nu¬ 
merically. However, near the bottom of the hyperbolic 
band where these shifts become small, they can be also 
found analytically. Thus, Eq. (25) gets replaced by 


qo - 


eo 


e-Lrf- 


2gtop + 29b ^ 


|ei| > eo,es 


(27) 


rlhot 


which shows explicitly that qo goes to zero as uj ap¬ 
proaches I where £■*■ sharply increases. 

Unlike in Fig. 4(a, d), in 4(b, c, e, f) the higher-order 
u > 1 modes are more difficult to see because of their 
lower relative intensity compared to those of the plasmon 
n = 0 (and n = 1) modes. Nevertheless, these modes 
remain well defined (underdamped). Near the bottoms 
of the respective hyperbolic bands their momenta qn still 
form an equidistant sequence except with a spacing 

27r 

q„+i - 9„ ~ , (28) 


which is modified compared to Eq. (23). This result can 
be obtained from Eq. (22) by approximating the finite 
differences such as </)top(qn+i) — 4’topiqn) by means of the 
derivative. Parameter I is defined by 


I = -2 


^4^top 

dq 


rf ^^bot 

dq 


(29) 


The physical meaning of this quantity is clarified in the 
next Section. 


IV. GOOS-HANCHEN EFFECT 

In this Section we consider the problem of the plasmon- 
polariton mixing from the point of view of real-space tra¬ 
jectories of the HP^ excitations. The question we con¬ 
sider is how polariton wavepackets propagate inside the 
slab and, in particular, how they reflect off its interfaces. 
As mentioned in Sec. I, for a given w, the angle 0 be¬ 
tween the z-axis and the group velocity v vector of HP^s 


is nearly independent of q. Therefore, monochromatic 
HP^ wavepackets propagate as highly directional rays. 
Naively, one would then expect that the polariton rays 
should zigzag inside the slab returning to each interface 
periodically with the repeat distance of 2d|tan0| = |(5|. 
Although such geometrical optics picture is adequate for 
insulating hyperbolic materials,®^ it is not quite correct 
for TI with gapless doped SS. The geometrical optics 
neglects a lateral shift or displacement of the rays af¬ 
ter each reflection [compare Figs. 5(a) and (b)], which 
is analogous to the GH effect of light. To define such 
a displacement one usually considers a wavepacket with 
a smooth envelope (for example, a Gaussian), in which 
case the displacement is the shift in the position of its 
maximum. 

While the GH effect^® was discovered measuring the 
reflection of light off an air-metal interface, the dis¬ 
placement 1 of the reflected ray is a general wave phe¬ 
nomenon^® that arises due to the dependence of the re¬ 
flection phase shift (j) on the lateral momentum q = 
{k^, k^). For example, the GH effect should also occur for 
surface plasmons.®^ The expression for 1 has the form®® 

1 = -SRe ^ . (30) 

oq 

It seems to be another general rule that the momentum 
dependence of (j) is significant only if the interface sup¬ 
ports electromagnetic modes with either a large propa- 



FIG. 5. (Color online) Polaritonic GH effect in TI slabs, (a) 
Schematics of the HP^ ray reflection in the absence of the SS. 
(b) The same with the SS. The wavy lines symbolize virtnal 
Dirac plasmons. The GH shift I is indicated, (c) The electric 
field distribution inside and/or at the upper surfaces of two 
slabs with equal <5 = —2.2d but different doping. The lower 
(“doped”) and the upper (“undoped”) parts of the image are 
computed for \p = a and 0, respectively. The split gate — 
a pair of metallic half-planes separated by a distance 2a — 
launches highly directional HP^ rays that bounce inside the 
slabs creating periodic “hot stripes” at their upper surfaces. 
The period is larger in the “doped” slab. The two small cir¬ 
cles, one in the undoped and one in doped part, are the rep¬ 
resentative locations of the HP^ reflections. Their enlarged 
views are shown in, respectively, (a) and (b). 
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gation length or a long decay length if such modes are 
evanescent. In the original photonic GH effect this is the 
case under the conditions of the total internal reflection. 
The magnitude |1| of the GH displacement can be inter¬ 
preted as the decay length of the evanescent transmitted 
wave. Alternatively, a large GH shift can occur if the 
interface supports surface plasmons or polaritons.®®^®® 
Experimental demonstration of the GH effect enhanced 
by surface plasmons of the air-metal interface has been 
reported.^® 

Comparing Eqs. (29) and (30), we recognize the length 
scale I in the former as the sum of the GH shifts due to 
the top and the bottom interfaces. Therefore, we con¬ 
clude that the Dirac plasmons must act as the transient 
interface modes for the HP^ rays bouncing inside the TI 
slab. Using Eqs. (20), (30), and taking into account that 
Ifte Cl ^ 9m ei, we find the GH shift at the top interface 
to be 


I 


top — 


4 

Qtop 


(eo 


9m Cl 

S)’ + hP' 


(31) 


A few comments on this result can be made. First, the 
GH shift is positive in our case, which means the dis¬ 
placement is in the same direction as the in-plane group 
velocity of the ray. Second, Ztop depends on the permit¬ 
tivity of the environment. For example, at fixed q, it 
vanishes if cq is very large. Conversely, for fixed eo, the 
GH shift reaches its maximum 


2 Apeo9mei ^ _ ^tt 

TP (Sfteei)^-I-(9mei)^ ’ ^ co^top ’ 


(32) 


at q = Tr/Xp. Finally, Zmax depends linearly on the char¬ 
acteristic size Xp of the Dirac plasmon wavelength and 
inversely on the absolute value |ei| » 9m ei of the effec¬ 
tive permittivity of the hyperbolic medium. 

In Fig. 6, we show Zmax for Bi 2 Se 3 and G/hBN/G sys¬ 
tems as a function of ui spanning their respective hyper¬ 
bolic bands. The relative shift, Imax/Xp, is greater in 
G/hBN/G because |ei| is smaller. Yet the absolute /max 
at the same = 0.3 eV is greater in Bi 2 Se 3 (where it is 
^ 200 nm) because it is hyperbolic at lower frequencies 
and Xp is larger at smaller to. 

One possible setup for experimental detection of the 
GH effect in TI is shown in Fig. 5(c). It differs from 
Fig. 1(b) in the addition of a split gate between the TI 
slab and the substrate. If this gate is made of a good 
conductor with large permittivity, it would suppress the 
GH shift at the bottom surface. However, it would serve 
another useful purpose. Previously, it has been demon¬ 
strated^^ that in the presence of an external oscillating 
field, thin metallic disks or stripes can launch HP^ in 
hBN. The split gate is to perform the same function here. 
The HP^s are preferentially emitted from the regions of 
highly concentrated field near the sharp metallic edges. 
We expect the rays to zigzag away from their launching 
points returning to the top surface with the period / — S, 
which is the sum of —6 « |5| due to the roundtrip in¬ 
side the slab and I = /top due to the GH shift at the 
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FIG. 6. (Color online) Maximum GH shift /max (in absolute 
units and as a fraction of Xp) for (a) TI slab and (b) G/hBN/G 
structure with the same chemical potential /r = 0.3 eV. 


top surface. Since / depends gtop, which is controlled by 
doping, the GH effect can be detected by measuring the 
positions of the electric field maxima [“hot stripes” in 
Fig. 5(c)] as a function of /r in the experiment. Although 
/ is quite small, the shifts accumulate after multiple re¬ 
flections, which can facilitate their detection, as in the 
original work of Goos and Hanchen.^^ 

To model the response of the system shown in Fig. 5(c) 
quantitatively we proceed as follows. We approximate 
the half-planes of the split gate by perfect conductors in 
the z = Q plane with the edges at a: = ±a. Let V(x, 0) be 
the scalar potential at z = 0 due to the external uniform 
field and all the charges induced on the gate. (Here and 
below the common factor is omitted.) Let V{k^) 

be the Fourier transform of U(a;,0). Using the notations 
for the reflection coefficients introduced in Sec. H, we 
express the potential V(x, z) inside the slab 0 < z < d 
by the integral 


V(x,z) 

t{kn 



V{k^)t{k^)e 


ik^x-\-i\k^ 1 2 tan 0 


1 - rio(fc^)ri2(fc^) 

1 - rio(A:“)ri2(A:“)e*'="'5 ’ 


(33) 

(34) 


For a consistency check we can consider the large-x be¬ 
havior of this inverse Fourier transform, which should be 
dictated by the poles of the integrand. These poles can 
be recognized as the HP^ momenta <?„ [Eq. (22)]. Since 
qn form the equidistant sequence [Eq. (28)], their super¬ 
position should indeed create beats of period 1 — 5, in 
agreement with our ray trajectories picture. Fig. 5(b). 

Explicit calculation of V{x, 0) requires a self-consistent 
solution of the Maxwell equations for our complicated 
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multilayer system, which is computationally intensive. 
Fortunately, very similar results for V(x, z) are obtained 
with little effort by approximating the true V{x, 0) with 
the “bare” potential that would exist in the TI is re¬ 
moved, that is, if c? = Ap = 0. At distances less than 
c/uj from the gap in the gate, this bare potential has the 
simple analytical form. 




+1; 

X < —a 

2 . / / A 


-arcsin(a:/a), 

\x\ < a, 

TT 


-1, 

X > a, 


(35) 


familiar from classical electrostatics. Its Fourier trans¬ 
form is given by 




(36) 


where J[j{x) is the Bessel function of the first kind and 
Vo is potential difference between the two parts of the 
gate. The tangential electric field corresponding to this 
potential. 


E 


X — 



(37) 


exhibits an inverse square-root divergence at the edges, 
which enables the localized HP^ emission. 

Carrying out the quadrature in Eq. (33) numerically, 
we have calculated the components and also the ampli¬ 
tude of the electric field E = yjE'^-\- E'^ over an interval 
of X a few 151 in length and z varying from 0 to d. Our 
results for E = E{x, z) for two doping levels, correspond¬ 
ing to Ap = 0 (undoped SS) and Xp = a (doped SS) are 
illustrated by the false color plots in Fig. 5(c). These 
plots are superimposed on perspective projections of the 
two slabs (doped and undoped), which are placed next 
to each other for easy comparison. The remaining pa¬ 
rameters of the calculations are S = —2.2d and a = O.ld. 
We see that a finite shift of the “hot stripes” at the top 
surface z = d exists in the doped case. This seems to 
vindicate our intuition but actually the situation is a bit 
more subtle. The problem is that the momentum distri¬ 
bution of our source [Eq. (36)] is very different from what 
we assumed it to be in the beginning of our discussion 
of the GH effect. This distribution is not narrow and 
not centered at some finite k^. Instead, it has positive 
and negative harmonics of equal strength and a long 
power-law tail at \k^\ ^ 1/a. The reason why the GH 
shift persists in our case is the spatial separation of the 
harmonics: due to the directionality of the HP^ prop¬ 
agation, the stripes to the left (right) of the launching 
points are created predominantly by negative (positive) 
kx- Since 1 has the same direction as q = (fc^,0), the 
stripes shift away from the origin on both sides of the 
y-axis. A formal derivation of this result can be done by 
splitting the integral in Eq. (33) into the k^ > 0 and the 
kx < 0 parts and identifying the relevant poles kx = Qn 
using contour integration methods. 


From numerical experiments with different a, we found 
that the largest shift of the stripes is obtained for a ^ Xp. 
This can be explained by arguing that the shift is maxi¬ 
mized when the characteristic ~ tt/o contributing to 
the integral in Eq. (33) is close to the momentum n/Xp 
at which I — in Eq. (31). 

Experimental detection of the “hot stripes” and their 
doping-dependent GH shift is possible via the s-SNOM 
imaging. This technique involves measuring the light 
scattered by the tip of an atomic force microscope 
brought to the sample and scanned along its surface. 
Using clever signal processing methods, it is possible to 
isolate the genuine near-held component of this scattered 
light, which originates from conversion of evanescent elec¬ 
tromagnetic waves emanating from the sample into free- 
space photons. In the proposed experiment, the evanes¬ 
cent waves are due to the HP^ modes launched by the 
split gate. The spatial resolution of the s-SNOM imaging 
is set by the tip curvature radius R. For typical R = 20- 
40 nm, it is barely sufficient to observe the predicted GH 
shifts in hBN/G, Fig. 6(b). Nevertheless, detecting the 
cumulative shift after several stripe periods should be fea¬ 
sible. The prior success of s-SNOM imaging experiments 
of surface plasmons and polaritons in graphene and hBN 
structures^°’^^“^®’^^’^®’^^’^® gives us a hrm conhdence in 
this approach. Note that if a doped graphene layer only 
partially covers the top surface of hBN, one literally gets 
the situation depicted in Fig. 5(c), where the doped and 
undoped regions are positioned side by side. 

In the case of Bi 2 Se 3 where the GH shift ~ 200 nm 
[Fig. 6(a)] is much larger, the spatial resolution of the s- 
SNOM is even less of an issue. The main obstacle is the 
scant availability of suitable THz sources. We are opti¬ 
mistic that in a near future this problem can be overcome 
as well. 


V. SUMMARY AND OUTLOOK 

Recent experiments®d4 have shown that coupling be¬ 
tween Dirac plasmons and bulk phonons of bismuth- 
based TIs should be strong. In this paper we have stud¬ 
ied this interaction taking into account the anisotropic 
phonon spectrum of such TIs. We have predicted that a 
TI slab can act as a tunable waveguide for phonon po¬ 
laritons, with the doping of the surface states being the 
tuning parameter. In additional to the change in disper¬ 
sion, the phonon-plasmon coupling can cause measurable 
real-space shifts of the polariton rays. Similar phenom¬ 
ena have been recently studied in artificial structures 
made by stacking graphene layers on top of hBN. The 
present work indicates that the TIs are a promising alter¬ 
native platform for realizing highly tunable, strongly con¬ 
fined, low-loss electromagnetic modes in a natural mate¬ 
rial. Additionally, while hBN/G waveguides operate in 
mid-infrared frequencies, Bi 2 Se 3 and similar compounds 
extend the same functionality to the technologically im¬ 
portant THz domain. 
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We envision several directions for further work in this 
field. One is to attempt a multi-source coherent con¬ 
trol of polariton emission and propagation using ultrafast 
laser pulses. A variety of such techniques has been devel- 
oped®° in the context of THz polaritonics of LiNbOs and 
LiTiOa. (Incidentally, a theoretical proposal®^ of inte¬ 
grating graphene into such materials would lead to polari¬ 
ton waveguides similar in functionality and perhaps also 
tunability to those studied in the present work.) Another 
intriguing direction is to explore oscillating spin currents 
which were predicted to accompany charge density cur¬ 
rents produced by Dirac plasmons.^^ It may be also in¬ 
teresting to study the effect of optical hyperbolicity^® on 
the high-energy bulk plasmons of the TIs.®^’®^ Finally, it 
may be worthwhile to investigate new applications that 
can be enabled by tunable hyperbolic polaritons. Har¬ 
nessing such types of modes for hyperlensing®^^®® or fo- 
cusing^^’^^ has been widely discussed. The present work 
shows that the GH effect and its dependence on doping 
and dielectric environment of the TI can be another av¬ 
enue for applications, for example, THz chemical sensing 
or characterization of spatially inhomogeneous TI sam¬ 
ples. We hope our work can stimulate these and other 
future studies. 
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Appendix A: Near-field spectra 

A fully realistic modeling of the s-SNOM imaging ex¬ 
periments proposed in Sec. IV is an unwieldy task re¬ 
quiring a repeated solution of the Maxwell equations for 
a system with complicated material properties, a hierar¬ 
chy of widely different length scales, and no special sym¬ 
metries. In this Appendix we present some results of less 
ambitious calculations that simulate a simpler structure 
depicted in Fig. 1(b). Although no split gate is present 
in this structure, the measured signal is still expected to 
reveal characteristics of the collective modes. In this case 
these modes are excited by the sharp tip itself. Hence, 
the tip plays the role of both the launcher and the detec¬ 
tor of the HP^ modes. Unfortunately, this implies that 
only the local response can be measured, which is a su¬ 
perposition of responses due to a distribution of momenta 
up to ~ 1 /i? rather than one specific q. 

We assume that the TI slab and the substrate are in¬ 
finite and uniform in x and y coordinates, so that the 
imaging capability of the s-SNOM is irrelevant. Instead, 
the quantity of interest is the frequency dependence of 
the measured near-held signal s(a;). A few more ex¬ 
planations about our calculational scheme are in order. 
We model the tip as a metallic spheroid with the curva¬ 
ture radius R = 40 nm and total length 720 nm. We use 


the quasistatic approximation but include the radiative 
corrections included perturbatively. This model®^’®® has 
been successful for simulating many recent s-SNOM ex¬ 
periments, and should be especially suitable in the THz 
domain where no antenna resonances or other strong re¬ 
tardation effects®® should appear. Our calculations in¬ 
corporate the so-called far-held factors, ®^^®® which are 
expressed in terms of rp(g, w) at g ~ w/c. This factors 
account for the fact that the incident wave is originally 
created by a far-held source and the scattered wave is ul¬ 
timately measured by a far-held detector. Finally, what 
we compute is not the full scattering amplitude s but its 
third harmonic S 3 , which is what experimentalists typi¬ 
cally report. The idea is that in the experiment the tip 
is made to oscillate at some low frequency D, so that s is 
periodic with this fundamental tapping frequency. The 
third Fourier harmonic of s, which is S3, gives a good 
representation of the genuine near-held signal. 

Naively, one can think of S 3 (a;) as a weighted average 
of the surface rehectivity rp{q,uj) over q. The weighting 
function has a broad maximum near q = qt, which in 
this case is equal to qt = 0.025nm“^ [the dashed lines in 
Fig. 4]. The presence of strong maxima of 3m rp due to 
collective modes with momenta q qt tends to enhance 
S 3 (a;). In a more rigorous picture,®® the maxima of 83 ( 10 ) 
correspond not to the resonances of the sample alone but 




a; (cm ^) 


FIG. 7. (Color online) Simulation of the s-SNOM signal S 3 
for BhSea slabs on a substrate with es = 10. (a) Fixed y = 
0.29 eV and different d. (b) Fixed d = 120 nm and different 

y. 
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to those of the coupled tip-sample system. The coupling 
can decrease the resonance frequencies by as much as®^^®® 
10-20 cm“^ compared to those seen in 9m rp maps. 

Our results for Bi 2 Se 3 slabs of various thickness d and 
chemical potential /r are shown in Fig. 7. Pairs of distinct 
peaks as well as smaller additional features are readily 
seen. In each trace, the stronger and sharper peak is lo¬ 
cated close to I = 64cm“^. The height of this peak 
decreases as d decreases [Fig. 7(a)]. However, its position 
is independent of d [Fig. 7(a)] or /x [Fig. 7(b)], which sug¬ 
gests that it is not related to the dispersive HP^ modes. 
Indeed, we have verihed that this prominent peak is al¬ 
most entirely due to the far-held factor |1 + rpp, which 
has a narrow maximum at ^ where rp « 1 . 

Each of the doped samples also produces smaller peaks 
in 53 ( 0 ;), of which the most prominent ones are those 
located near uj = 146cm“^ and w = 163cm“^, the upper 
boundaries of regions B and C of Fig. 3. The position and 
especially the strength of the peaks is /x-dependent. As /x 
increases, the peaks grow in height and gradually shift to 
higher frequencies, see Fig. 7(b). These peaks are due to 
the surface modes: the n = 1 mode of region B and the 
n = 0 mode just above region C, see Fig. 4(d, e). The 


increase of the peak heights with fi can be qualitatively 
explained by the increase of the absolute value of rp. 
The shift in position is unfortunately more difficult to 
interpret without a better understanding of the effective 
weighting function that relates 9 mrp(( 7 , ox) to 53 ( 0 ;). 

While /X > 0 traces are due to combined action of plas- 
mons and phonon-polaritons, the /x = 0 one is expected 
to reveal the phonon-polariton response. Interestingly, 
that trace exhibits a sharp dip at w = 163 cm“^, see 
Fig. 7(b). We have checked that this dip is not caused 
by the far-held factor. However, its relation to the HP^ 
modes of Fig. 4(d) is not obvious to us. 

The thickness dependence of S 3 is illustrated in 
Fig. 7(a). As one can see, the near-held peak at a; = 
163 cm“^ has a broad high-frequency side, which system¬ 
atically expands as d decreases. This trend rehects the 
blue shift of the n = 0 mode dispersion in thinner slabs, 
compare Fig. 4(e) and (f). 

Overall, our simulations predict that the near-held re¬ 
sponse of Bi 2 Se 3 slabs should exhibit systematic spectral 
changes with doping and thickness that are measurable 
by the s-SNOM. Such experiments may provide insights 
into properties of tunable HP^ modes of these novel sys¬ 
tems. 
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